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We propose and illustrate an approach to coarse-graining the dynamics of evolving net- 
works (networks whose connectivity changes dynamically). The approach is based on the 
equation-free framework: short bursts of detailed network evolution simulations are coupled 
with lifting and restriction operators that translate between actual network realizations and 
their (appropriately chosen) coarse observables. This framework is used here to acceler- 
ate temporal simulations (through coarse projective integration), and to implement coarse- 
grained fixed point algorithms (through matrix- free Newton-Krylov GMRES). The approach 
is illustrated through a simple network evolution example, for which analytical approxima- 
tions to the coarse-grained dynamics can be independently obtained, so as to validate the 
computational results. The scope and applicability of the approach, as well as the issue of 
selection of good coarse observables are discussed. 
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INTRODUCTION 



Complex dynamic systems, exhibiting emergent dynamics, often arise in the form of graphs (or 



networks): the internet, social networ 



tion networks and more 



IE 



iS, chemical and biochemical reaction networks, communica- 
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30-32 



36l | . In a social network, for example, the individuals 



are represented by nodes (or vertices) , while the relations among them are represented by the edges 
connecting these nodes. 

One type of network dynamics arises in cases where the network topology (connectivity) is 
static, but the state of each vertex is a variable that evolves in time (in part due to interactions 
with the states of connected vertices) . Such problems are said to exhibit "dynamics on networks" . 
A different type of network dynamics (on which we focus here) arises when the existence or the 
strength of connections between the vertices constitute the variables that evolve in time. These 
problems are said to exhibit "dynamics of networks". These two types of dynamics are not, 
of course, mutually exclusive; clearly we can have dynamical problems involving dynamics both 
on and of networks. The evolutionary network problems we will refer to in this work will be 
exclusively of the second type of dynamics mentioned above: dynamics of networks, where the 
network structure is the state that changes over time. We will restrict ourselves to networks with 
unweighted edges, that are either present or absent; we will not study edges of continuously variable 
strength here, even though our methods can be adapted to function in such situations also (in fact, 
with appropriate modifications, in any type of network evolution problem). 

In our networks of choice the detailed graph representation constantly changes over time accord- 
ing to some specified rules: edges (and/or nodes) are added or deleted. Although these fine-scale, 
"node level" , microscopic details of the graph are changing in time, coherence often emerges at a 
macroscopic level. At such a coarse-grained, system scale, the (expected) structure is often seen to 
evolve smoothly in time: it may eventually become stationary or may possibly oscillate between 
a number of states. We will use a coarse-grained, macroscopic view to study graph dynamics, 
treating the coarse temporal evolution as a dynamical system. 

Reduced descriptions of high-dimensional dynamical systems (in our case, of large graphs) are 
possible predominantly due to two mechanisms: (a) decoupling of the evolution of a set of variables 
and/or (b) separation of time scales between evolution of different groups of variables. In the 
former (much simpler) case, the few uncoupled variables evolve independently of other variations in 
the system and hence a reduced closed description can be written in terms of just these variables. 
In the latter (more interesting) case, characteristic time scales of evolution of a few variables 



3 



(called the "slow" variables) of the system are much longer than the time scales of evolution of the 
remaining "fast" variables. After a short evolution time the fast variables will then often become 
slaved to the slow variables, i.e., the evolution of the fast variables becomes (in the long term) 
solely determined by the evolution of the slow variables. The long-term dynamics of the system 
can therefore in principle be approximated by equations written only in terms of the slow variables, 
the "coarse variables" of the system. Note also that, depending on the time scales of interest, it 
may be possible to close the system (write closed form evolution equations) at different levels of 
detail. 

To use the established tools and techniques of dynamical systems (e.g. bifurcation and stability 
computations), however, one must explicitly derive the evolution equations describing the coarse 
variable dynamics. We are interested in cases where such coarse evolution equations conceptually 
exist, yet they are not available in closed form. In such cases, it has been shown that it may 
be possible to circumvent their explicit derivation using equation- free techniques 31. These 
techniques are based on short bursts of appropriately designed computational experiments with 
the detailed, "fine scale" network dynamic evolution model, and on the knowledge of the appro- 
priate coarse observables: the variables in terms of which reduced closed coarse equations could in 
principle be written. Using traditional numerical algorithms as the basis for the design of com- 
putational experiments, and exploiting algorithms that translate between coarse variable values 
(relevant network statistics) and actual realizations of networks with these statistics, equation-free 
approaches may significantly accelerate the computer-assisted study of network dynamics. In re- 
cent work, we have applied equation-free techniques to illustrative examples of several different 

r n n 

types: molecular dynamics [10[, collective animal behavior [28^, cell population dynamics [7] and 



also dynamical models on static networks Here, we demonstrate the use of equation- free tech- 
niques on an illustrative graph evolution problem, and test our results against explicitly derived 
coarse equations. 

A crucial prerequisite for equation-free modeling is the knowledge of a good set of coarse vari- 
ables - the collective network features that can be used to predict its (expected) coarse evolution, 
the variables in terms of which the coarse model would close. While a large graph is an intrinsi- 
cally high dimensional object and difficult to visualize, its complicated structure can be probed by 
measuring statistical properties of the graph. Such statistical properties of graphs are often good 
candidate coarse variables. Commonly used statistical properties for describing a graph include the 
average degree, the degree distribution, the clustering coefficient, and the distributions of shortest 
path lengths. [30,] . 
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As we will see below, even after a good set of such coarse variables has been selected, an 
important requirement for using our equation-free algorithms is the ability to routinely construct 
graphs exhibiting prescribed values of these properties. It is clearly trivial to construct a graph 
with a given number of nodes and edges; yet it is quite non-trivial to construct realizations of 
graphs with, say, a prescribed distribution of shortest path lengths. 

These issues will be discussed and illustrated through a simple model example in the remainder 
of the paper as follows: The model behavior is discussed briefly in Section|Tll We describe our coarse- 
graining approach in Section IIIII For our simple example, it so happens that several theoretical 
results can be explicitly obtained (and used for validation of the computer-assisted approach). 
These will be discussed in Sections IIVI and |Vl We will conclude with a brief summary and a 
discussion of the scope of the approach, its strengths and shortcomings, and of certain (in our 
opinion) important open research directions. 

II. MODEL: A RANDOM EVOLUTION OF NETWORKS 

We consider a simple, illustrative model of stochastic network evolution. Let the graph at any 
discrete time, T = 0, 1, ... be denoted by G„(T). The subscript n denotes the number of nodes in 
the network. The rules governing the dynamics at each iteration can be described as follows: 

1. Select a pair of nodes at random and connect them together by an edge if they are not 
already connected to each other. 

2. With probability r, remove an edge chosen uniformly at random, (r stands for removal 
probability.) 

We performed numerical simulations using these detailed, node-level, "microscopic" rules (using 
r = 0.9) on graphs with n = 100 nodes and observed the evolution of several statistical graph 
properties over time. In these preliminary numerical experiments, the initial conditions were either 
empty graphs or Erdos-Renyi random graphs with a specified value of edge probability, p (of 
which empty graphs are a special case, corresponding to p = 0). It is interesting to consider the 
evolution of graph properties starting from an ensemble of initial conditions: Erdos-Renyi graphs 
with the same edge probability p. Fig(T] shows the evolution of one such property of interest: the 
degree distribution of the (nodes of the) graphs, obtained statistically from a 100-copy sample. 
(The degree of a node in a network is the number of edges connected to the node.) The figure 
demonstrates that the degree distribution can be thought of as a function (ignoring, for the sake of 
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FIG. 1. The (apparently) smooth evolution of degree distribution according to the model with r = 0.9. 



simplicity, the discrete nature of the degree) that appears to evolve smoothly over time. Thus, the 
degree distribution can serve as a potential coarse observable of the system; one must, however, 
carefully investigate whether it is a good candidate for coarse variable. 

In particular, one must test whether it is possible to obtain a description of the long-term 
evolution of this observable that is closed - that would imply that an accurate reduced model of 
the system evolution can, at least in principle, be derived. In other words, if one had measurements 
of the chosen set of coarse variables (observables) at particular time, that information alone should, 
in principle, be enough to predict future states of the system (in terms of these variables). We 
reiterate that, depending on the time scales of interest, different useful reduced models of the same 
system, i.e. closures at different levels of coarse description, may be possible to obtain. 

We then proceeded to test whether the degree distribution constitutes a good choice of coarse 
observable. For this purpose we constructed initial graphs with identical degree distributions but 
different higher order information (triangle statistics, for instance), evolved the graphs using our 
model rules and compared their evolutions in time. Figures [2)^a) and [2]|^b) show the evolution of 
degree distributions and of clustering coefficient distributions starting from graphs with identical 
degree distributions, but distinct distributions of clustering coefficients. (The clustering coefficient 
of a node is the ratio of the number of triangles attached to the node divided by the maximum 
possible number, given its degree.) The solid (blue) curves represent the case where the initial graph 
is an Erdos-Renyi graph (p = 0.25). The dotted (red) curves corresponds to the case where the 
initial graphs were created to match the degree distribution of Erdos-Renyi graphs with p = 0.25; 
this was done by sampling n = 100 numbers (degrees) from the required degree distribution and 



using the Havel-Hakimi 



161 ] algorithm to construct a graph with the sampled sequence of degrees. 



The Havel-Hakimi algorithm consists of three iterated steps: 
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FIG. 2. Evolution of (a) degree distribution and (b) triangle distribution, for two distinct ensembles of 
graphs. The solid (blue) curves represent the case where the initial graphs are chosen to be Erdos-Renyi 
graphs with p — 0.25. The plots corresponding to this case are denoted by the labels, El, E2 and E3 at 
times 0, 10 and 20 respectively. Note that one time unit in the plot corresponds to n = 100 iterations of the 
model. The dotted (red) curves represent the case where the initial graphs were created to match the degree 
distribution of the initial graphs of the previous Erdos-Renyi case through the Havel-Hakimi algorithm. The 
plots corresponding to this second case are denoted by the labels, HI, H2 and H3 at times 0, 10 and 20 
respectively. 

1. sort the vertices by their degrees (dj) in non-increasing order; 

2. select the first vertex (di) and connect it to the next di vertices; and 

3. decrease the value of di by di (it is now 0) and the value of the di successive degrees by 1. 

As an illustration, consider the sorted list of residual degrees of the vertices (residual degree is the 
degree of a vertex minus the number of assigned edges for the same vertex at any given point in 
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the algorithm) 

di > d2 > ■■■dn- 

Once we connect the first vertex with the next di vertices, the residual degree sequence is 

(0, d2 - 1, - 1, ddi+i - 1, c?di+2, •■•)• 

After sorting, the sequence becomes 

{d2 - 1, ds - 1, dd^+i - 1, 0). 

These steps are repeated until we get a sequence of zeros, implying that the graph with the desired 
degree sequence has been achieved. 

For both our test cases, we thus have the same initial degree distribution. However, since the 
graphs themselves are different, properties like clustering coefficients will differ between them, at 
least initially. The model was evolved using 100 copies to obtain smooth distribution functions. 
The figures show that the evolution of degree distribution is similar in both cases. Although 
the initial clustering coefficient distributions are very different for the two cases, they quickly 
(within a few time steps, compared to the time scale of evolution of the degree distribution) 
converge visually to the same function and subsequently coevolve. This suggests that the clustering 
coefficient distribution (and possibly other higher order information) may become quickly slaved 
to the evolution of the degree distribution: triangle statistics (and so also clustering coefficients) 
evolve at a much faster rate, and quickly reach a distribution that appears to depend only on the 
comparatively slowly evolving degree distribution. These results encourage us to attempt to find 
a coarse-grained reduction of the system using a discretization of the degree distribution as the 
coarse variable(s). 



III. COARSE-GRAINING 



We propose a computer-assisted coarse-graining approach -the Equation-Free (EF) framework 
Q) Q]~ to develop and implement a reduced description of the system, using the degree distribu- 
tion as the coarse observable. In this approach, short bursts of simulations at the "microscopic" 
(individual node) level are used to estimate information (such as time-derivatives) pertaining to 
the coarse variables. This is accomplished by defining operators that allow us to translate between 
coarse observables and consistent detailed, fine realizations. The transformation from coarse to 
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fine variables is called the lifting operator {R), while the reverse is called the restriction operator. 
If we denote the microscopic evolution operator by (pt, the macroscopic evolution operator can be 
defined as 



As an illustration, we implemented coarse projective integration |19l . |20| ] using the equation- 
free approach and the degree distribution as the coarse variable. In coarse projective integration, 
the system is integrated forward in time using occasional short bursts of detailed, microscopic 
simulations at the level of individual nodes, and the results are used to estimate time derivatives 
at the level of the degree distribution. In the following discussion, time units are taken to comprise 
n = 100 iterations of the model. We used Erdos-Renyi random graphs {p = 0.25) as the initial 
conditions and repeatedly performed the following operations: 

1. Simulation: The model is executed initially for 10 time steps (in each time step, we perform 
n = 100 iterations of the rules of the model). 

2. Restriction: The graph degree distributions, fi{d), are observed from simulations period- 
ically (at intervals of 2 time steps). We stored the degree distribution at discrete values 
/.ii = /i((ii), where di = 0, 1, 2, ...99. The distribution may also be discretized using a smaller 
number of points and interpolated at intermediate values (we will discuss the possibility of 
alternative, more parsimonious representations below). In fact, our plots of degree distribu- 
tions are smooth interpolations from this 100 discrete value representation. 

3. Projection forward in time: The time-series of the coarse variables over the final segment 
of the simulation (in our case, the last 3 observed discretized degree distributions) are used 
to predict the new values of the coarse variables at a future time (in our case, 10 time steps 
down the line). This is done on the basis of established numerical integration schemes (in 
our case, simple forward Euler). For the simple model differential equation dx/dt = f{x), 
the forward Euler scheme would read 



The difference in our case is that the time derivative (/'(x(t)) above) does not come form a 
closed formula, but is instead estimated from the recorded recent time series segments. In 
our computations, this time-derivative estimation and projection in time is performed as fol- 
lows: Given the last 3 observed discretized degree distributions, the associated 3 cumulative 




x{t + At) w x{t) + Atf'{x{t)). 
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distribution functions are found, and the degrees, gi, corresponding to uniformly distributed 
percentile points, pi G {0, 0.01, ...1}, i = 1, 101, such that 

rgi 

/ fi{x)dx = pi. 
Jo 

Thus, the pairs of points {gi,Pi) constitute discrete approximations of the cumulative degree 
distributions. The values of these percentile degrees, gi, observed at the last 3 time steps, 
are the variables we actually projected forward in time, estimating the time derivative of 
the corresponding forward Euler scheme by fitting a straight line, and extrapolating it for 
10 further time steps. When projecting the discretized version of cumulative distributions, 
one should take care to retain monotonicity of the predicted (projected) distributions [39|. 



In our simulations we did not encounter such a numerical problem, possibly because (a) we 
used several copies to get smoothened degree distributions, and (b) the projection times were 
relatively short. 

4. Lifting: To restart the simulation, the predicted discretized distribution must be trans- 
formed into consistent graph realizations. We accomplish this by using the Havel-Hakimi 
algorithm; we may have to sample the projected degree distribution until we draw a graphi- 
cal sequence of degrees. (A graphical sequence is a sequence on non- negative integers that is 
realizable as a degree sequence of a graph.) Checking if a sequence is graphical is performed 
as a part of the Havel-Hakimi algorithm: if the algorithm terminates successfully, the se- 
quence is graphical; otherwise, it is not. In the latter case another sequence is sampled until 
a graphical one is found. Once we have these graphs, the procedure repeats: we continue 
simulations for 10 more time steps as in stage 1, keep the last 3 observations, and so on. 

The results of coarse-projective integration are shown as (blue) thin lines in Fig. [3] (note that 
degree distributions are plotted only when we collect simulation data, and not when we jump in 
time). The results from direct integration of the full model are plotted (red, thick lines) every 
20 time steps for comparison; the evolution of the degree distribution predicted by our coarse 
projective integration clearly coincides with the full detailed simulation. 

We studied the rate of (temporal) convergence of discretized degree distribution as a given 
sample network evolution approaches steady (stationary) state. Fig. U] shows the evolution of the 
difference between the instantaneous degree distribution fj,{d) close to the steady state and the 
steady state degree distribution /i°^(d) itself. We used PC A to evaluate the first few principal 
components of the time sequence (//(d) — n°°{d)). The first two singular values were found to be 
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FIG. 3. Coarse projective integration (CPI) with the degree distribution as the coarse variable. Resuhs 
from CPI are plotted as (blue) thin lines every 2 time steps during simulations, while results from direct 
simulation are plotted as (red) thick lines every 20 time steps for comparison. Note that one time unit in 
the plot corresponds to n = 100 iterations of the model. 

~ 1.00 and 0.0175 respectively. The vectors corresponding to these first two singular values are 
plotted as (blue) dots in Fig. [5j These eigenvectors represent directions along which the decay 
of the degree distribution, from the given randomly chosen initial condition to the steady state, 
is the slowest. We found that the principal components were qualitatively similar for a wide 
variety of choices of initial graphs. For the specific example shown in the figure, the two principal 
components were found to be well-approximated by the analytical expressions e^^~^^^^ ^"^^{x — ll)/5 
and e^^~^^^'^ {x — 8){x — 14)/20, which are plotted as (red) solid lines in the same figure. The 
relevance of the qualitative form of these principal components will be discussed later in Section |Vl 




FIG. 4. Decay of the degree distribution close to steady state: evolution of the difference between the 
instantaneous degree distribution /^(d) and the steady state degree distribution ^°°{d). 



11 




FIG. 5. The first two PCA components of the ensemble of vectors — fi°°{d) shown in Fig. |4] are plotted 
as (blue) dotted lines. These components were found to be well-approximated by the analytical expressions 
e(^-ii)^/20(2, _ ll)/5 (first component, left) and e'^^~^^^^^'^^{x — 8){x — 14)/20 (second component, right), 
plotted as (red) solid lines for reference. 

In addition to coarse projective integration, we also performed coarse fixed point computa- 
tions. Instead of finding the (discretized) stationary degree distribution (coarse fixed point of the 
evolution) through direct simulation, one can also obtain it by solving the equation 

F{fi{d)) := fi{d) - <^Md)) = 0, 

where is the coarse time-stepper over t time steps. We find the roots of F using a damped 
Newton-Krylov GMRES iteration scheme Is^. The standard Newton algorithm updates the 
value of ^(d) by /i„ 4.1(d) = /i„(d)+A^((i), where Afi{d) satisfying [DF{fin{d))]Afi{d) = —F{fin{d)). 
Since ^tipid)) is not explicitly available (but can be evaluated through the coarse time-stepper), 
its Jacobian cannot be obtained by analytical differentiation; in the Krylov-based approach the 
action of this Jacobian on known vectors (its matrix-vector products) are estimated through nu- 
merical directional derivatives. Thus, linear (and through them, nonlinear) equations are solved 
through iterative matrix-free computations; these methods are naturally suitable for equation- free 
computation, where explicit Jacobians are not available in closed form. 

There are a couple of points worth mentioning about the use of this general methodology in 
the case where the unknowns solved for constitute a (discretized) distribution function n{d): (a) 
the iJ,{d) vectors that arise should be non-negative and (b) they should integrate to 1. At every 
iteration of the root-finding algorithm, these two properties should be preserved. These conditions 
on fi{d) can be stated as conditions on A^(d): Afj,{d) < —fin{d) and Afi{d) should integrate to 
0. We satisfy the first condition by using a damped Newton-Krylov method. The correction term 
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FIG. 6. Convergence of the Newton-GMRES iterations in steady state computations: the norm of the 
function to be solved for via Newton-GMRES is plotted against the iteration number. 



A^(d) is scaled by a constant to ensure non-negativity of the result: fin+iid) = fin{d) + clS.^[d) 
for c sufficiently small. The second condition is guaranteed by the structure of the Krylov method 
itself. When solving [Z)F(/i„((i))]A/i((i) = —F{fj.n{d)) using a Krylov method, the solution Aii{d) 
lies in the space spanned by • [Z)F(^„(d))]" • for 

some finite m < n. Because each of the spanning vectors has integral 0, the update term A/Li((i) 
(and hence cAfi{d), for any constant c) also has integral 0. The convergence of this procedure can 
be shown by plotting \\n{d) — (j){fi{d))\\ at every iteration as in Fig. |6l 

Coarse projective integration and coarse-fixed point algorithms are only two illustrations of 
equation-free computation: even though explicit coarse-grained evolution equations are not avail- 
able, the assumption that they in principle exist helps solve the coarse-grained model through 
appropriately designed short bursts of detailed simulation (also through lifting and restriction). 
Many additional computational tasks (e.g. coarse continuation and bifurcation computations, 
coarse eigencomputations, coarse controller design and even optimization) also become possible in 
this framework [2l||. Our computations so far provide numerical corroboration of the possibility of 
coarse-graining our network evolution model: a reduced model appears to close (accurately enough 
to be usefully predictive) in terms of only the (discretized) degree distribution. In what follows, we 
will provide certain theoretical results to support our choice of coarse variables and provide some 
intuition about the overall coarse-graining approach. 



IV. THEORETICAL JUSTIFICATION 



In this Section we discuss theoretical results that corroborate the observed separation of time 
scales between the evolution of our coarse variables (the discretized degree distribution) and higher 
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order information about the evolving graph statistics (e.g. the distribution of triangles). Such 
information would support our ability to usefully close a coarse-grained model at the level of 
degree distributions. 

For this simple graph evolution model, theoretical results for the evolution of edge density and 
vertex degrees can be easily derived using basic notions of probability; this was one of the reasons 
for choosing this model as our illustration. In general, for obtaining such results (including results 
for the evolution of triangles and more), one makes use of the notion of dense graph limits, which 
will be outlined later. 

Recall that our model graphs evolve in discrete time steps according to the rules given in Section 
im Gn{T) denotes the graph in n nodes at any discrete iteration, T = 0, 1, ... Gn{T, is the entry 
in the corresponding adjacency matrix representing the edge between nodes i and j. Let E{GniT)) 
represent the set of edges in the graph. 



A. Evolution of edge density 

We denote by e„(T) = \E{GniT))\ the number of edges in the graph at a given discrete time, 
T. Let Pn{T) be the edge density of the graph, GniT) at time T: 

Note that the evolution of en{T) is itself a Markov chain, and that it is decoupled from the other 
variables: 

P(e„(r+l)-e„(r) = l|G„(r)) = (l -p„(T)) • (l - r), (1) 
P( e„(T+l) - e„(r) = I Gn{T) ) = (1 - Pn{T)) ■ r + p„(T) • (1 - r), (2) 

p(e„(r+i)-e„(r) = -i|G„(r)) =p„(r))T. (3) 

In order to study the evolution of Pn{T), we introduce the continuous time variable t £ [0,oo) 
and let T = [Qt\ , so that T + 1 corresponds tot + dt with dt = ■ Let p„(t) := p„([(2)tj)). 
It fohows from (dJ, ([2]), ([3]) that we have 

E( p„(t + dt) - pn{t) I pn{t) ) = (1 - r - pn{t))dt, (4) 
T>'^{pn{t + dt)- pn{t)\pn{t)) = {{l- Pn{t))Pn{t)+r-{l-r))(^\ dt. (5) 
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FIG. 7. Coarse projective integration (CPI) of the edge density evolution: Results obtained from CPI using 
edge density as the coarse variable are plotted as blue dots (see text for details). The trajectory obtained 
from full direct simulation (red solid line) is plotted for reference. One time unit in this plot corresponds to 
n — 100 iterations of the model. 

Letting our process evolve for T >i n'^ steps (i.e. t x 1), we obtain 

E(pn(t) -/9„(o)) ir-dtxi, 

D^(p„(i)-p.(0))iT.Q^'dtxn-l 

Since the variance of Pn{t) — Pn{0) in the above formulas is much smaller than its expected value 
if 1 ^ n, we can replace Pn{t) by E(/9n(t)) in (|4]) without causing significant error. This leads to 
the deterministic equation 

Ap(t) = (1 - r) - Pit) (6) 
for p(t), the limit Pnit) p{t) as n — )■ oo, corresponding to ([4]). Thus 

p(t) = (l-r) + (p(0)-(l-r))e-* (7) 

and limt^oo pit) = l-r. 

From dH, it is clear that the evolution of edge density p{t) is truly decoupled from the evolution 
of other quantities - one does not need the higher order statistics of the graph to get slaved to it 
over a short time period on a sort of slow manifold for the reduced model to close. As mentioned 
earlier, this represents one type of conditions for which reduced descriptions of the system become 
possible. In Section IIIII we computationally implemented a reduced model using the entire degree 
distribution as the coarse observable(s). ([6]) represents an even more reduced description that 
meaningfully closes in terms of edge density. Similar to the degree distribution case, we implement 
coarse-projective integration using now only edge density as the coarse variable. The single scalar 
value of the edge density was observed and recorded during brief bursts of direct simulation; the 
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time horizons for direct simulation and subsequent projection forward in time were again chosen 
to be 10 time steps each. The hfting operation in this case involves creating graphs with a given 
value of edge density. Our particular implementation of this created graphs where every edge was 
selected to be present with a probability equal to the edge density. The results of coarse projective 
integration with the empty graph as initial condition, shown as (blue) dots in Fig. [71 agrees visually 
with the evolution of edge density as observed from direct simulations, shown as a (red) solid line. 
This corroborates the closure of a reduced description in terms of edge density only. We now 
proceed to derive results supporting our previous observation of closure of a reduced description 
at the more detailed (less coarse) level of degree distribution. 



B. Evolution of the degree of a vertex 

Consider the time evolution of the degree of a node i £ [n], i.e. the number of other vertices i 
is connected to. The order of magnitude of the degree of i is n. Define a normed degree, Dn{t), of 
a vertex i at scaled time t as 




we omit the z-dependence from the Dn{t) notation for the sake of simplicity. Following a derivation 
similar to the one used in Section [IV A|, we obtain: 



E(l?„(t + dt) - Dn{t) I Gn{t) ) = [l-(l + Dn{t)^ dt, (9) 

{Dnit + dt) - Dn{t) I Gn{t)) = U - Dn{t) + r^^) fl-'dt. (10) 

Now (jlOp is of smaller order than ([9|) if n S> 1, so that we may replace Dn{t) by E(Z)„(t)). By 
analogy with the argument of subsection IIV Al we see that Dn{t) — >• D{t) as n — >• oo, where D{t) 
solves the ODE 

Substituting the explicit formula ([7|) into (jlip and using the formula for the solution of inho- 
mogenous linear ODEs, we obtain an explicit solution for D{t): 

D{t) = p{t) + eT^p{t)^ (^p(0)^L>(0) - /9(0)i^) . (12) 

Clearly, lim^^oo D{t) = limt^oo p{t) = 1 - r. 
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Note that the evolution of the degree of a given vertex depends both on its current degree and 
the current edge density of the graph as a whole. Yet, if the degree distribution of a graph is given, 
the edge density of the graph also follows from it. This clearly supports our observation that the 
system closes at the level described in Section Hill in terms of the degree distribution. Approximate 
differential equations were derived so far to describe the evolution of (expected values of) the 
normed degree and of the edge density of graphs. The next step is to explore the dynamics and 
influence of higher order information, like triangles. In order to derive similar results for statistics 
of triangles, we take advantage of concepts from the theory of convergent dense graph sequences, 
of which a rigorous and formal introduction can be found in 



C. Convergent dense graph sequences 

Let Gn{T) denote the adjacency matrices of our evolving graphs, and let the individual entry 
representing the existence of an edge between nodes i € [n] and j G [n] at time T be denoted by 
Gn(T, The notion of dense graph limits will be useful for describing the time evolution of the 
statistical properties of Gn{T) when n ^> 1. The limit object of a convergent graph sequence is a 
so-called graphon W, where W : [0, 1]^ — )■ [0, 1] is a measurable (but not necessarily continuous) 
function with the properties W{x,y) = W{y,x) and W{x,x) = 0. Assume that for each n G N we 
have a graph G„ with vertex set [n]. We now informally define the notion of convergence of the 
sequence Gn to W, i.e. G„ — )• W. 

One might heuristically imagine the adjacency matrix of Gn as a black-and-white television 
screen (a white pixel at position represents an edge between vertices i and j); a convergent 
graph sequence becomes then a sequence of TV sets showing the same picture at higher and higher 
resolution. The limiting graphon W will then be the picture seen on the "perfect TV" where each 
point {x,y) G [0,1]^ is a "pixel of infinitesimal size"; the local density of black/white pixels will 
then give us the impression of shades of grey. For the precise definition of the so-called cut-distance 



25l . Theorem 



S{Gn, W) between a finite graph and a graphon, see [25|, l26|. Section 4]. Note that by 
6.13] every graphon W arises as a limit for a convergent graph sequence G„. 

Clearly, there exist many adjacency matrices corresponding to different labelings of the nodes 
of the same graph, and in the definition of Gn W we are allowed to relabel the vertices of Gn 
(i.e. to rearrange the pixels our TV set). Correspondingly, we are allowed to relabel [0,1] using 



measure-preserving transformations in order to obtain equivalent forms of the graphon W (see 



Section 3.1]). For the purposes of the present paper, accounting for rearrangements is not required. 
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If F is the adjacency matrix of a small graph on k nodes, then we define the homomorphism 
density t{F,Gn) by 

t{F,Gn):=j^_ E 2[Vz,iG[A;]:F(i,i)<G„(v^(i),V^(j))], (13) 

where we sum over all possible injective functions (/? from [k] to [n]. F is our test graph and t{F, Gn) 
the homomorphism density of F in G„. 

We define the homomorphism density of in by 

t{F,W) = [\ . . [' n W{x,,xjf^'^^^dxi...dxn. (14) 
■^0 -^0 l<i<j<k 

Denote by Kk the complete graph on k vertices; for example, K2 is an edge and is a triangle. 
Erasing an edge from a triangle gives a "cherry", a simple graph with three vertices and two edges. 

\F(C^ )\ 1 

Now, denote by pn := Ls.' the edge density of Gn- It follows from [25!, Section 6.2] that 

[2) 

Gn^W implies t{F, Gn) ^ t{F, W): 

pw ■■= lim pn= lim t{K2,Gn)=tiK2,W)= [ [ Wix,y)dxdy, (15) 

hm t{Ks,Hn) = t{Ks, W) = C C [' W{x, y)W{y, z)W{z, x)dxdydz, (16) 

Jo Jo Jo 



lim t(cherry, if„) = t(cherry, 1^) = / / / W{x,y)W{y, z)dxdydz. (17) 

Jo Jo Jo 

It is the ability to write such equations that makes working with graphons useful for our purposes. 
Once the graphon is identified, one can approximate the density of any test graph F in Gn, 1 <^ n 
using expressions similar to Equations p^. p^ . pT|) . 



D. Evolution of the graphon 

If we consider a convergent graph sequence Hn — )• W, where Hn is a graph on n vertices, and 
for each n we run our Markov process Gn(T) with initial state G„(0) = Hn, then 

Gn{[('l)t\)^Wt, 



2^ 

where Wt{x,y) is the solution of the following ODE: 



^Wt{x,y) = l- (^l + -^]Wt{x,y). (18) 



The heuristic derivation (1201) of this formula is based on 



Wt{x,y)^P{Gn{[(^fjt\,[x-n\,[y.n\) = l) (19) 
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where T = [Q)*] , dt = Q) ^ , i = [x ■ n\, j = [y ■ n\. Note that e(T) p{t) (2) . 



Wt+dt{x,y) - Wt{x,y) = 

P{G{T + l,i,j) = l\G{T,i,j)=0)-{l-Wt{x,y)y 
P{G{T + l,i,j)=0\G{T,i,j) = 1)-Wt{x,y) = 

2] ■i^-Wtix,y))-^yWtix,y) = 



l-[l + ^^jWt{x,y))dt. (20) 
This results in (jlSp . 

By substituting the exphcit formula ([7]) into (jlSp and using the integral formula for the solution 
of inhomogenous linear ODEs, we obtain an explicit solution for Wt{x,y): 

Wt{x,y) = p{t) + eT^pit)^ (^p{0)^Wo{x,y) - p{0)^^ . (21) 

We have limt^oo Wt{x,y) = p{t) = 1 — r. Thus for t S> 1 the graphon becomes almost 

constant with value 1 — r; (jl9p shows that the stationary state of our graph dynamics looks like 
an Erdos-Renyi graph with edge density 1 — r. In addition, (jlSp implies that an initially constant 
graphon will remain constant at future times. Thus, the family of Erdos-Renyi graphs is an 
invariant set of the dynamics of system. 

Note that, given the explicit formula for Wt{x,y) we may obtain an explicit formula for the 
density of triangles and cherries in Gn (T) using (jl6p and (jl7p . Differential equations for describing 
the evolution of triangles and cherry densities at time t can be found directly by differentiating 
those equations. 

E. Evolution of cherry and triangle densities 

Recalling (fT7|) we have 

i (cherry, Wt) = [ [ [ Wt{x,y)Wt{y,z)dxdydz. 
Jo Jo Jo 

Differentiating (fT7|) with respect to time and using the graphon evolution result of p^ . we get 



d_ 
di 



2p{t) - 2 ( 1 + ) t(cherry, Wt). (22) 
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An ODE describing the time evolution of the density of triangles (in terms of the edge density and 
the cherry density) can be similarly derived: 

^t(i^3, Wt) = 3 t(cherry, Wt) t{K^, Wt). (23) 

F. Convergence rates 

Now that we have equations describing the statistics of degrees, cherries and triangles, we can 
find out the rates at which these quantities converge to their steady (stationary) states. A function 
a{t) converges to another function b{t) at rate a if 

lim - log(|o(t) - b{t)\) = -a. 

t— >CJO t 

In order to establish this, it is enough to prove that 

t->oo \a[t) — b[t)\ at 

From d?]) we can directly see that the edge densities of our graphs converge to the steady state 
value of 1 — r at a rate a = 1. This can also be shown by using ([6]) to write the following equation: 

±(p{t)-{l-r)) = - (pit) - {1 - r)) . 

We now show that the normed degree D{t) of a vertex converges to p{t) at a faster rate. From 
(©,(111]) we obtain 



A (D{t) - Pit)) = -(^1 + ^^ (Dit) - Pit)) . (24) 

Thus a = lim(_).oo (^1 + — this case. Note that a = also follows from the explicit 

formula (jl2p . For example, if r = 0.9 then a = 10. 

From (fT8|) . we similarly obtain that for any x,y & [0, 1] the function Wtix,y) converges to pit) at 
a rate a = 

If the graphon Wt evolves according to (fTS|) and pit) = Wtix,y)dxdy then pit) solves ([U]). 

If we let Wtix,y) = pit) for any x,y £ [0,1] then Wt also solves (fT8|) . Thus, the set of constant 
graphons is invariant under the dynamics. We now show evidence that this "invariant manifold" 
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is actually attracting: 




(25) 



This implies that t(cherry, Wf) converges to t(cherry, Wf) at rate a 



2j^. If r = 0.9 then a = 20 



for this case. Thus, if we evolve graphs that already possess the steady state values of their edge 
density, their cherries will converge to their steady state value twice as fast as the rate at which 
normed degrees evolve to their corresponding steady state. 

We now consider triangle density evolution. Let W'^ and W'^ be two distinct graphons with the 
same values of edge and cherry densities: 

t{K2,W^) = t{K2, and t(cherry, H^^) = t (cherry, 

If we let and evolve according to the dynamics (jlSp with initial states and W"^, 
respectively, then by ([6]) and ([22]) we have 



for all t > 0. The fact that the density of edges and cherries coincide for and "helps" the 
densities of triangles in and Wf to converge to each other even more rapidly: 



Thus the rate of convergence of the relative triangle density is a = 3^;^: for r = 0.9 this works out 
to be 30. This result, in particular, explains why we observed an apparent slaving of the triangles, 
as discussed in in Fig. [2l The (blue) solid and (red) dotted curves there showed the evolution of 
two graphs with the same degree distribution. Graphs with the same degree sequence also have 
the same number of edges and cherries, which implies (j26p . The number of triangles in these two 
graphs (corresponding to the two cases in Figl2]) converge to each other three times as fast as the 
rate at which the degree distributions ultimately evolve. 

We now estimate convergence rates from direct simulation results, confirming the validity of our 
approach and approximations even for relatively small systems: in these results n = 100 nodes and 
we simulate the model using a value of 0.9 for the parameter r. Figure [8] shows how information 
about degrees, cherries and triangles converge to their corresponding steady (stationary) states. 



t{K2, Wt') = t{K2, Wf), t(cherry, PF/) =t(cherry, W'^ 



(26) 



dt 




21 




Time.t Time.t Time.t 



FIG. 8. Evolution of the logarithms of quantities related to (a) degree, (b) cherries and (c) triangles. The 
initial conditions for the three cases are explained in the text. SND denotes the sum of squares of the normed 
sequence of degrees of all the 100 nodes. Nch denotes the number of cherries in the graph. Nt,er and 
Nt.hh denote the number of triangles in two graph evolutions starting with an Erdos-Renyi random graph 
with p = 0.1 and a graph created using the Havel-Hakimi algorithm respectively, both with the same degree 
distribution. Note that one time unit in this plot corresponds to (2) = 4950 iterations of the model. Slopes 
of the three lines are —10.44, —18.49 and —29.98 respectively, which are in close agreement with theoretical 
results of —10, —20 and —30 respectively. 

Note that, in all these cases, the time t is scaled so that one time unit comprises (2) iterations of 
the model. Logarithms of quantities (defined in the caption of the figure) related to these statistics 
are plotted in the y-axis versus time. For producing the first two plots (corresponding to degrees 
and cherries), the initial graphs were created by first sampling from a normal degree distribution 
with mean 10 and standard deviation of 1. For comparison, the steady state graphs have a degree 
distribution whose mean and standard deviation were evaluated to be ~ 10 and ~ 3 respectively. 
Thus, the initial graphs have the same mean degree (and hence the same edge density) as the steady 
state graphs, but differ from these steady state graphs in their actual detailed degree distribution. 
Since the graphs are initialized with the steady state edge density, this edge density remains close 
to its steady state value at all times. From ()24|) and ()25p . the convergence rates for the quantities 
in the first two figures are expected to be 10 and 20 respectively. This successfully matches the 
numerical convergence rates (the absolute values of the slopes of evolutions) of the terms related 
to degrees and cherries in the figure: they are seen to be 10.44 and 18.49 respectively. For the 
third plot, containing information about triangles, we simulate the model from two different initial 
graphs with the same degree distribution but different number of triangles. The first simulation is 
initialized with an Erdos-Renyi random graph with p = 0.1. For the second case, we created initial 
graphs using the Havel-Hakimi algorithm, using the degree sequence of the first case as input. 
Since graphs with the same degree sequence also have the same number of edges and cherries. 
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(j26p is satisfied. Hence from (j27p . we expect a convergence rate of 30 for the relative number of 
triangles between these two graphs, which successfuhy matches the numerically computed value of 
29.98. Thus, although the theoretical results are in principle accurate only at the limit of very large 
graphs, all the numerical values we computed using graphs with only 100 nodes are remarkably 
close to the limiting theoretical values. 



V. SOME ADDITIONAL THEORETICAL RESULTS 
A. An SDE for the degree of a vertex 

In the previous section we approximated the evolution of the normed degree through a deter- 
ministic ODE, arguing for the relative smallness of the order of magnitude of the variance of what 
is really a stochastic process. In order to now describe the evolution of the stochastic process Dn{t) 
at a finer level, we approximate it by a stochastic differential equation (SDE) rather than an ODE: 



dD^it) =(l- Dn{t) - r^) dt + Jl- Dn{t) + • n-^'^dWu (28) 

V Pn\t) J y Pn\t) 

where Wt now denotes the standard Brownian motion. To rationalize the choice of such an approx- 
imation, we observe that the solution of (|28p satisfies ([9]) and (jlOp and that the Dn{t) defined by 
([8|) satisfies \Dn{t + dt) — Dn{t)\ < ^, i.e. the trajectory of Dn{t) is very close to being continuous. 
One can then suggest that it is appropriate (and even natural) to use the Brownian motion in ()28p 
as a driving function. For rigorous results validating the SDE approximation of discrete stochastic 
processes, see jl?! ]. 

We are interested in the fluctuations of Dn{t) around its expected value Dn{t) = E(Z)(t)). 
Using ([9]) we can see that Dn{t) approximately solves the ODE 

dDnit) = ^1 - Dn{t) - ry^^ dt. (29) 

If we define X{t) := n^/'^{D{t) - Dn{t)) then by ([MD and ^ X{t) solves the SDE 



dX{t) = -(l + X{t)dt + Dn{t) + r^^dWt. (30) 

If n is big enough, then Pn{t) ~ p(t) and Dn{t) ~ Dn{t) ~ D{t), so we may use the deterministic 
p{t) and D{t) in the right-hand side SDE of X{t) without causing much error. 

Since p{t) and D{t) are known and explicit (c.f. ([7]), (fT2]) ). (pOl) is a linear SDE, i.e. an Ornstein- 
Uhlenbeck process with time-dependent drift and diffusion coefficient. From this, it follows that 



23 



X{t) can be approximated by a Gaussian process with mean and an explicit formula for the 
variance at time t. If we let t — )• c« then p{t) — )• 1 — r and D{t) — t- 1 — r, so (I30p becomes 

dX{t) = -j^X{t)dt + V2^dWt, (31) 

an OU process. The variance of the stationary distribution of this Markov process can be expressed 
using the drift and diffusion coefficients and it is normally distributed: 



j =J\f (m = 0,a = a/(1 - r)A 



X(oo) ~AA |^m = 0,fT2 = ^-^^ ) =AA(m = 0,f7 = v^(r^^7)7) (32) 

Now Dn{t) = Dn{t) + n-^/'^X{t), from which we get: 

D„(oo) ^M(m = l-r, a = -^\/r{\ - r)) (33) 



It is worth noting that similar results can be derived for the edge density by defining Yify := 
(2)^^^/0(0 - ^(*) solves the SDE 



dy(t) = -Y(t)dt + ^{\- p(t))p{t)^r-(\-r)mt (34) 

This is analogous to ([30|) . If we let t — >• 00 then p{t) — >• 1 — r, so becomes the following OU 
process: 



dy(t) = -y(t)dt + V2r • (1 - r)mt. 



B. A PDE for the evolution of the normed degree probability distribution 

If we consider the SDE ([3T]) and denote the probability density function of X(t) by P(t, x), then 
P(t, x) solves the Fokker-Planck (or Kolmogorov-forward) equation: 

|p(t,x) = ^A(^.p(t,^)) + ,^P(t,^). (35) 

A simple argument is that, when 1 ^ n, the trajectories of the evolving degrees of vertices i and j 
show little correlation, since the source of randomness for the degree evolution for distinct vertices 
is almost disjoint: they have at most one edge in common. It then follows that observations of the 
time evolution of the empirical degree distribution histograms (for the entire graph) may be well 
approximated by solutions of the Fokker-Planck equation ([35|) (for the degree probability density 
of a single node). 
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Solving the eigenvalue-eigenfunction problem corresponding to the differential operator on the 
right-hand side of ()35p gives rise to the Hermite differential equation, whose eigenfunctions are the 
Hermite functions. In particular, the first three eigenvalues and eigenfunctions are 

Ao = 0, /o(2;) = exp (^-i;::^!^) 

A2 = -TTF, /2(x)=exp(-i^:^) -(r- j^) 

• foix) is (a constant times) the density function of ([32]) which is also the stationary solution 
of ([35]). 

• fi{x) represents the direction along which the decay of a generic initial distribution P(0,x) 
to P(oo,x) is the slowest: 

^lim exp(^— — ) • {P{t, x) — P{oo, x)) = c - /i(x) 

• /2(x) represents the direction along which the decay of P{0,x) to P{oo,x) is the slowest, if 
P(0, x) is a generic even function of x. 

These formulas corroborate the plots in Fig. [5l even though the assumptions made in order to 
derive the theoretical results are valid only at the limit of very large graphs. It is interesting to note 
that the leading principal components of the decay of the empirical degree distribution histograms 
to steady state, which we found to be well approximated by the functions e^^~^^'''^ ^^^{x — ll)/5 and 
g(a;-ii)2/20^^ — 8){x — 14)/20, Can also be very well matched to the Fokker-Planck eigenfunctions 
and /2(x) by shifting and rescaling coordinates. 

VI. SUMMARY AND CONCLUSIONS 

In this paper, we have demonstrated a computational framework for coarse-graining evolution- 
ary problems involving networks. We illustrated our methodology using a specific model with 
simple, random evolution rules. The proposed methodology applies, in principle, to any network 
evolutionary model with an inherent separation of time scales between the evolutions of a few 
important coarse variables, and the remaining slaved variables (observables) . For our illustrative 
model, we were able to analytically derive certain theoretical results, justifying our choice of coarse 
variables and quantifying the observed time scale separation. We used the notion of dense graph 
limits to formulate some of our arguments for successful computational coarse-graining. It is con- 
ceivable that some of the theoretical tools used here might be useful in deriving insights in other 
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dynamic network models. We emphasize, however, that for the right problems our coarse-graining 
procedure will work irrespective of whether one is able to analytically derive such supporting the- 
oretical results. 

The generality of the approach raises other important general questions in the area of complex 
networks. We mentioned earlier that a critical step is the identification of suitable coarse variables. 
There are at least two open questions that need to be addressed in that regard: 

1. How does one find the appropriate coarse variables for a given model? 

2. Once suitable coarse variables are identified, how does one solve the problem of constructing 
networks with prescribed values of the chosen variables? 

The second question is more concrete, and hence, we will address it first. For certain specific 
properties, there exist well-known algorithms to construct graphs with specified values of those 
properties. For instance, for constructing networks with a given degree distribution, we repeatedly 



used the Havel-Hakimi algorithm 
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161 ]. Alternative approaches to construct graphs with a 



model 



3a], the Chung-Lu 



specified degree distribution include, among others, the configuration model 

11], and a sequential importance sampling algorithm Q]. Beyond the degree distribution, 



there are only a few properties for which standard approaches have been established for constructing 



graphs with specified values of those properties. For example, algorithms that create grap 



the following properties can be found in the literatur e: g iven degree-degree distribution 
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IS with 
, given 



degree distribution and average clustering coefficient 221] , and given degree distribution and degree 
dependent clustering [35]. For other properties (or combinations of properties), more generalized 
mathematical programming approaches (e.g. [14]) could potentially be used to solve the network 
generation problem. 

The first question, however, requires more new ideas. For the example presented here, coarse- 
graining was originally based on computational model observations, and preceded the derivation 
of theoretical results. In general, the motivation for good coarse variables can come from standard 
heuristics, intuition about the model under consideration or observations of evolution of statistical 
quantities of the model. However, smart use of data mining tools such as diffusion maps 



d require the definition of a useful metric 
9|, Is^] ) . Automatic data mining to extract 



might provide answers in a more generic sense. This wou^ 
quantifying the distance between nearby graphs (see e.g. 
good coarse variables from model observations is, in some sense, a holy grail of model reduction 
methods. 
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